HUTlONAl.  BUREAU  OE  STANOAROS 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


ts. . 

lO 

(£)  . 


D  D  C' 

:  '  I — 1  nn  n  n  r 


f'r?  :',^r?nnnE| 

<  HOV  1  5  1379 

llLLiL^lialJ  U  LLl 

A 


SICuAlTv  CLAttl^iCATlOM  0*  TMlt  »*QC  0«« 

1  REFOirr  OOCUMCNTATION  PACE 

READ  tNSTRUrnnNS  1 

■trORC  COMPLKTTNU  rORM  | 

1  mumIII - 

4  T|Tl,|  (mMKU) 

The  Use  of  the  Empirical  Probability 
Generating  Function  to  Estimate  the 

Neyman  Type  A  Distribution  Parameters 

T  ftuTMoara, 

Harold  R.  Bishop 

•  .  COmTAaCT  OM  MuilSCf^ai 

t  •ca»ea«i«a  oao^hiiatiom  h*«c  4mo  AOoatM 

Naval  Postgraduate  School 
•Monterey,  California  93940 

•>  COMTaOkLIMO  0'*lCC  aAMC  4H0  400atM 

Naval  Postgraduate  School 
.Monterey,  California  93940 

OC^ItT  OATC 

September  1979 

«UM«CA  OR  »*OCt 

86 

'4  MOaiTeaiHO  4aiMCv  MAmI  4  400ac*tf»  rnllmfmmi  Imm  CannallMj  Olllfl 

Naval  Postgraduate  Sc.hool 
•Monterey,  California  93940 

•tCuOi'^V  CLAfti  f*  •*ttm 

Unclassified 

>4  ai*Tai»uTioa  tT4Tt>HMT  !■»(  mit  Htmari) 

Approved  for  public  release;  distribution  unlimited 

'T  Oi(Tai*uTieii  IT»Ttm«T  '•»  *•  tttirmi  aiMvW  »•  m—»  f0.  II  rnilmmmi 


>•  (u»aLCliCMr«aT  nortt 


't  BtY  aOMO*  rCmimmi  tm  MM*  il  «»nn»»  m»»  imrniilt  ^  mm*  aiMaari 


M  aMTMaCT  (CmttHmm  tm  fwna  II  minia  itmntr  ar  Mm*  «a*an 

The  Maxinium  Likelihood  estimators  for  the  N’eyman  Type  A  dis¬ 
tribution  parameters  are  very  difficult  to  compute.  In  this 
thesis,  the  Empirical  Probability  Generating  Function  is  used 
to  provide  estimators  that  are  easier  to  compute  and  have 
asymptotic  efficiency  at  least  as  high  as  97%  of  that  for  the 
Maxi.mum  Likelihood  estimators  over  most  of  the  parameter  space 
considered.  The  estimators  found  by  this  method  are  consistently 


»< 

I  JAM  rt 

(Page  1) 


00 


1473 


tMTiOM  0*  *  DOV  ••  i«  OMOklTI 


MCuWTV  CkAMI 


»iC<lTiow  0»  TwH  iW>i  ftf  >iw«| 


Approved  for  public  release; 
distribution  unlimited 


THE  USE  OF  THE  EMPIRICAL  PROBABILITY 
GENERATING  FUNCTION  TO  ESTIMATE  THE 
NEYMAN  TYPE  A  DISTRIBUTION  PARAMETERS 

by 

Harold  Ralph  Bishop 
Lieutenant,  United  States  Navy 
3.  A.,  San  Jose  State  College,  1971 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 

MASTER  OF  SCIENCE  IN  OPERATIONS  RESEARCH 

from  the 

NAVAL  POSTGRADUATE  SCHOOL 
September  19''9 


Author 

Approved  by: 


2 


ABSTRACT 


The  Maximum  Likelihood  estimators  for  the  Neyman  Type  A 
distribution  parameters  are  very  difficult  to  compute.  In 
t.his  thesis,  the  Empirical  Probability  Generating  Function 
IS  used  to  provide  estimators  that  are  easier  to  compute  and 
have  asymptotic  efficiency  at  least  as  high  as  97%  of  that 
for  the  Maximum  Likelihood  estimators  over  most  of  the  param¬ 
eter  space  considered.  The  estimators  found  by  this  method 
are  consistently  better  than  the  Method  of  Moments  and  the 
Method  of  Zero  Frequency  estimators  with  respect  to  asymptotic 
efficiency.  The  considerations  of  preference  in  using  one 
method  over  anot.her  are  discussed. 
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I. 


INTRODUCTION 


Several  'contagious'  distributions,  including  the  Neyman 
Type  A  distribution,  were  derived  by  Neyman  [1]  in  the  course 
of  modeling  the  reproduction  behavior  of  certain  types  of 
plants  and  bacteria.  More  currently,  these  models  come  under 
the  heading  of  branching  processes.  Several  methods  of  pa¬ 
rameter  estimation  have  been  proposed  for  the  Neyman  Type  A 
distribution  with  varying  degrees  of  success.  T.he  two  most 
frequently  used  are  the  Method  of  Moments  and  t.he  Method  of 
Zero  Frequency.  The  performance  of  an  estimation  scheme  can 
be  measured  by  the  asymptotic  efficiency  of  its  estimators 
as  defined  by  Wilks  [2].  This  paper  proposes  a  met.hod  for 
estimating  the  parameters  of  the  Neyman  Type  A  distribution 
by  using  the  Empirical  Probability  Generating  Function. 

For  t.he  Neyman  Type  A  distribution,  t.he  Method  of  Moments 
provides  estimates  that  are  quick  and  easy  to  compute,  al¬ 
though  the  efficiency  of  t.hese  estimators  is  not  always  high. 
It  is  well  known  that  Maximu.m  Likelihood  estimators  are 
asymptotically  efficient,  but  they  are  very  difficult  to  com¬ 
pute  for  the  Neyman  Type  A  distribution. 

The  method  of  estimation  by  using  the  Empirical  Probabil¬ 
ity  Generating  Function  provides  estimators  with  asymptotic 
efficiency  close  to  that  of  the  Maximum  Likelihood  estimators 
over  most  of  the  parameter  space  considered.  Extensive  tables 
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are  presented  along  with  an  algorithm  that  is  relatively  easy 
to  use  on  current  programmable  calculators. 

The  estimators  found  by  applying  this  method  have  con¬ 
sistently  higher  asymptotic  efficiency  than  the  Method  of 
Moments  and  the  Method  of  Zero  Frequency  estimators.  However, 
larger  sample  sizes  may  be  required  to  cause  the  algorithm  to 
converge  on  a  unique  pair  of  estimates. 

A  brief  outline  of  the  report  is  presented  as  follows.  A 
derivation  of  the  Neyman  Type  A  distribution  with  a  list  of 
some  of  Its  mathematical  properties  is  contained  in  section 
II.  The  equations  for  the  Maximum  Likelihood  estimators  and 
asymptotic  efficiency  appear  in  section  III.  The  mathemat¬ 
ical  details  of  using  the  Em.pirical  Probability  Generating 
Function  to  provide  estimators  are  found  in  section  IV,  along 
with  a  pictorial  display  of  asymptotic  efficiency  for  these 
estimators.  An  algorit.hm  for  computing  the  estim.ates  of  this 
met.hod  is  presented  in  section  V.  The  method  of  using  the 
Empirical  Probability  Generating  Function  is  compared  to  the 
Method  of  Moments  in  Section  VI.  The  algorithm  described  in 
section  V  IS  illustrated  by  a  n-omerical  example  along  wit.h 
computer  studies  in  section  VII.  Tables  for  use  with  the 
above  algorithm  are  found  in  section  VIII  and  a  summary  of 
the  results  appears  in  Section  IX. 


NEYMAN  TYPE  A  DISTRIBUTION 


II . 


Neyman's  Type  A  distribution  can  be  used  to  model  decay 
and  reproduction  processes  which  fit  the  assumptions  of  com¬ 
pound  Poisson  processes.  A  derivation  of  this  distribution 
provides  an  understandi.ng  of  its  structure,  usefulness,  and 
properties.  For  these  reasons  it  is  presented  i.n  the  fol¬ 
lowing  form. 

Consider  a  decay  process  in  which  radioactive  nuclei 
decay  to  produce  daughter  nuclei  according  to  a  Poisson  proc¬ 
ess  with  rate  '  .  Also  suppose  that  the  n'unber  of  daughter 
nuclei  produced  after  eac.h  decay  are  i.ndependent  ra.ndom 
variables  havi.ng  a  common  Poisson  distribution  wit.h  parameter 
mT .  Let  N(t)  =  the  n’umber  of  decays  in  the  time  interval 

(0,t] . 

Y,  *  the  number  of  daughters  produced  from  the  ith 
decay . 

X  =  the  total  number  of  first  generation  daughters 
produced  in  the  ti.me  interval  (0,t]. 


Then 


X 


N(t) 

i=l 


Because  N(t)  is  a  Poisson  process  with  rate 

-\t  (\t)” 

?  [N(t)  »  n]  *  e  n!  for  n  *  0,1, 


*9 


Conditioning  on  N(t)  is  used  now  to  find  the  probability  mass 
function  of  the  random  variable  X.  For  fixed  N(t),  X  is  a 
Poisson  random  v'ariable  with  parameter  nm2  when  N(t)  =  n. 

Thus 


*  nit  2  ( n  It  2 ) 


?  {X  =  X  1  N(t) 

=  n)  =  e 

x! 

for  X  = 

1,2,3, . 

The  event  that  no 

daughters  are 

produced 

in  the 

first  gener- 

ation  may  be  desc 

ribed  by  two  cases,  namely 

1)  X  =  C  If 

N  ( t)  =  0  or  , 

2)  X  =  0  if 

N  ( t )  ^  0  . 

From  case  1) 

?  ( X  =  0  1  N  ( t ) 

=  0)  = 

1  IS  obvious. 

- 

nm- 

And  from  case  2) 

?(X  =  0  ,  N(t) 

II 

o 

)  =  e 

Now,  using  the  total  probability 

theorem 

f 

p  =  ?(x  =  0) 

o 

=  :  PCX  =  0 

n=0 

N  ( t )  = 

n)  ?(N( 

t)  =  n) 

or 

-  t  » 

-nn2 

e  e 

-■•t  (Xt) 

n 

o 

II 

e  Z 

n=  1 

n ! 

-m^ 

-'■t(l  -  e  ‘) 

o 

II 

e 

• 

It  is  also  obvious  that  for  X  * 

X  >  0,  PCX  =  X  1 

N (t)  *  0)  * 

And  by  using  the  same  conditioning  argument  as  above  it  is 


evident  that 


=  P(X  =  X) 


ao 

r 

n=l 


xP 


-Xt  (Xt)“ 
e  n ! 


These  two  formulas  may  be  combined  into  one  expression  which 
IS  given  in  Shenton  [3]  as 

P  =  e  ^  m,^  [  0^  -t-  =1^  +  +  .  ]  (2.1) 

_ 2_  -yt-  -yr- 

where 

-m_  "  1  if  X  =  0 

2  X  I 

m.  =  • t  3  =  m.e  and  O'  =  < 

^  ^  0  if  X  *  0 

w 

This  probability  distribution  was  originally  derived  by 
Neyman  and  is  called  the  Neyman  Type  A  distribution.  Other 
forms  oi  equation  (2.1)  may  be  found  in  Johnson  and  Kotz  [4]. 

Some  of  the  properties  of  this  distribution  are  needed  in 
the  further  development  of  this  paper  and  are 
1.  Probability  generating  fu.nction. 


X 


E(t  )  =  e 

2.  Moment  generating  function. 


m-  (t-1) 

(  1  -  e  ] 


(2.2) 


uX 

E(e  ) 


■m^  [  1  -  e 


m^  (e  -  1) 


(2.3) 


By  using  the  well-)cnown  formula 


.  u.\ 

d*^  E(e  ) 
du-'^ 


0 


E(X'' 


1 


he  following  are  obtained 


(2.4) 


3. 

E(X) 

= 

4. 

E{X^) 

^  C  X 

X  ^ 

+  +  mj^m.,) 

5  . 

E(X^) 

=  ( 1 

X  ^ 

+  3m2  ■*■  3mj^m2  *  3nij^m2^  ■*■ 

m2^ 

2  2, 

^2  ) 

6.  VAR(X)  =  (1 


(2.5) 

(2.6) 

(2.7) 
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EFFICIENCY  AND  MAXIMUM  LIKELIHOOD  EQUATIONS 


A  method  for  estimating  the  parameters  of  a  distribution 
should  yield  estimators  with  high  asymptotic  efficiency  and 
It  IS  well  known  that  Maximum  Likelihood  estimators  have 
asymptotic  efficiency  equal  to  1.  However,  the  Maximum  Like¬ 
lihood  equations  may  be  difficult  to  solve  and  such  is  the 
case  for  the  Neyman  Type  A  distribution.  Replacing  som.e  of 
these  equations  in  the  Maxi.mum  Likelihood  system  may  produce 
a  system  of  estimating  equations  that  is  more  readily  solved 
and  still  produces  estim*ators  with  high  asymptotic  efficiency. 
Estim.ating  equations  considered  here  are  those  that  equate  a 
statistic  with  Its  expected  value.  This  paper  presents  one 
such  system  for  the  Neym.an  Type  A  distribution. 

It  can  be  shown  (see,  for  'xample.  Read  [5])  that  the 
asymptotic  efficiency  of  the  esti.mators  of  suc.h  a  substitute 
system  is 

.  1 

rpF  =  _ 

where  j  .1  =  information  determinant  of  the  probability  dis¬ 

tribution  considered  and 

T 

M  =  lim  InE(e-eM'‘-d)  ] 
n 

with  •?  “  vector  of  distribution  parameters 
and  d  =  vector  of  esti.mators  of  -. 
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e  above  expression  for  EFF  can  be  simplified  if 


h  =  vector  of  estimating  equations 


A  =  the  matrix  E  ^  ~ 

'  **  *  rs  -♦  »x> 


nE  (hh"^) 


EFF  = 


(3 


!■••  c 

The  Maxim.um  Likelihood  equations  for  the  Seyran  Type  A 
stribution  are  derived  in  [3]  and  the  follcwinu  form.s  are 
)und  in  i -i  I 

1)  X  =  m,m^ 

r,  3 


2)  -  (X.  -  1)  -^1  =  nx 


13. 


’.ere  x  =  t.he 

sam.p  le 

mean 

nd  n  =  the 

sam.p  le 

sice. 

It  should 

be  not 

ed  tha 

licated  rune 

tion  of 

•"’l 

f  some  other 

set  of 

equat 

'.s  as  substitutes  for  the  equa- 
lons  in  n.2),  it  is  necessary  to  com.pute  the  inform.aticn 
eterminant  of  t.he  Seyman  Type  A  distribution.  The  infcrma 
ion  determinant  has  been  derived  in  [3]  and  there  it  is 


hewn  that 


IV 


ESTIMATING  EQUATIONS 


The  Maxi.Tiun  Likelihood  equations  are  quite  difficult  to 
solve  for  the  Neyriuin  T%’pe  A  distribution  as  seen  from  the 
equations  in  (3.2).  Namely,  the  second  equation  is  the  one 
that  produces  the  complication.  Replacing  the  second  equation 
in  (3.2)  with  a  different  one  will  allow  for  simpler  calcu¬ 
lations  if  t.his  new  equation  is  chosen  judiciously.  However, 
the  idea  of  maintaining  high  efficiency  should  also  be  con¬ 
sidered  when  making  this  c.hoice  of  a  substitute.  The  use  of 
the  Empirical  Probability  Generati.ng  Fu.ncticn  (hereafter 
referred  to  as  EPG"'  see.m.s  to  be  appropriate  in  attaining 
this  goal.  This  was  suggested  by  the  successful  application 
of  this  choice  in  the  case  of  the  Negative  Bincmial  distri¬ 
bution  as  demonstrated  by  Caglayan  [T].  Since  the  EPGF 
contains  the  indepe.nde.nt  variable  t,  the  solution  of  t.he  new 
system  of  estimating  equations  can  be  "tuned"  so  as  to  ac.hieve 
t.he  esti.mators  of  ma.xi.m.um.  efficie.ncy  capable  from  t.he  new 
system.  Considering  the  form  of  the  probability  generating 
fu.nction  for  the  Neym.an  T-ype  A  distribution,  it  is  easily 
amendable  to  programming  on  ourrent  programmable  calculators. 
Wit.h  this  m.ctivation  in  mind,  the  first  Maximum  Likelihood 
equation  in  (3.2)  and  the  EPGF  equation  which  follows  are 
considered  as  a  possible  alternative  to  estimiating  t.he  pa¬ 
rameters  for  t.he  Neyman  Tj'pe  A  distribution. 
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The  EPGF  is  defined  as 

n  X 

EPGF  =  1  :  t  ^ 

n  1*1 

where  n  *  the  sample  size  and  =  the  ith  sample  value.  The 
expected  value  of  the  EPGF  is  the  probability  generating  func¬ 
tion  of  the  distribution  of  the  random  variable  X  and  is  given 
by  equation  (2.2).  The  system  of  estimating  equations  con¬ 
sidered  above  is  then 

1)  X  =  m.m^  (4.1) 

1  2 

m^  (t-1) 

n  X.  -m,  [  1  -  e  ] 

2)  1  :  t  ‘  *  e  ^  (4.2) 

n  1  =  1 

The  notation  m,  *  and  m^*  is  used  to  denote  t.he  estimators 
1. 

of  m,  and  m.,  found  by  solving  equations  (4.1)  and  (4.2).  To 
solve  these  equations  the  variable  t  must  have  a  value  in  the 
interval  [0,1)  to  insure  that  t.he  sign  of  the  left-hand  side 
of  equation  (4.2)  is  compatible  with  the  right-hand  side. 

The  value  c.hcsen  for  t  (whic.h  will  be  denoted  as  t*)  is  found 
bv  ma.ximizing  the  asymptotic  efficiency  of  t.he  esti.matcrs  m,  * 
and  m.,*. 

To  com.pute  the  asymptotic  efficiency  of  m,  *  and  m^* 
equation  (3.1)  is  used.  To  reduce  the  amount  of  typing  in 
the  following  expressions  let 

m^  (t-1) 

-m,  I  1  -  e  ] 

G(t)  =  e  ^ 
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The  vector  of  equations  (h)  is  given  by 

h. 


X  ni  ^  2 


The  matrix  A  is 


*=  1 

n 

1 

X 

n 

1  =  1 

3h^ 

3h^ 

3m2 

£ 

3  hi 

3h2 

3m^ 

3  m2 

that 

Its 

el  erne  n 

11 


-m. 


-m. 


1  ^ 


m^ ( t- 1 ) 


G(t) 


(t-1) 


A-,  T  =  m,  ( 1-t)  e  G  ( t)  . 


The  matrix  C  is 


lim.  nE 
n 


f  vj  ^  h  h  ] 


:  h^h,^  h^‘ 


and  It  follows  that  its  elements  are 
^11 


lim  nVAR(x) 
n 


mj^m2(l  *  ^2) 


^12 


n  X. 

lim  nCOV(x,  1  It) 


COV(  x^,  t  ) 


(4.3) 

(4.4) 

(4.5) 

(4.6) 


(4.7) 


n— n  i»l 
the  last  step  resulting  from  x^  and  x  being  independent  when 


1  X  :. 
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Continuing , 


A  • 

Ci2  =  E(x^t  *■)  -  E(x^)E(t  M 

and  using  the  fact  that  for  well-behaved  distributions  (in¬ 
cluding  Neyman's  Type  A) 

X  X  . 

E(x  t  =  tE  (  _d(t  1  =  tC  (t) 

dt 


then 


(t-1) 

^12  *  nij^n2G(t)  I  te  *■  -  1  ].  (4.8) 

3y  symr.etry  C,  ^  =»  C,,. 

**  1  ^  X 


1*: 

^  1 

n  X 

X 

C-5-,  =  lim  n\ 

'A.=i(i  :  t  ^) 

=  VAR(t  ^) 

(4.9) 

**  **  p» 

n  1  =  1 

thus 

C22  =  ^(7‘) 

-  lG(t)l^. 

(4.10) 

Using  equations  (4.3)  -  (4.6) 

m.,  (t-1) 

A,  =  :n^G(t)  [  1  *  (m^Ct-l)  -  l)e  ^  1  (4.11) 

and  using  equations  (4.7)  -  (4.10)  with 

^^(t-l) 

3  =  1  m.,  nij^ni-,(te  “  -  I)*" 

then 


Cj  ■  [  ( l+m, )  G  ( t^ )  -  G“(t)  B  ]. 

i,  tL  ^ 


Now  substituting  (4.11)  and  (4.12)  into  (3.1)  results 

2  m2(t-l)  2 

m.G  (t)  (  If  (m-t  -  -  1)  e  ^  'i 

ZFF  »  — - - - i - 

|Alm^  (  (l*m^)G(t‘)  -  G"(t)  B  ] 

^  lb 


(4.12) 


in 


(4.13) 


1' 


t*  IS  the  value  of  t  used  to  solve  equations  (4.1)  and  (4.2). 
It  is  defined  as 

EFF(t*)  =*  max  EFF(t). 

0<t<l 

In  order  to  use  this  definition  the  values  of  and  m2 
must  be  given.  t*  and  EFF(t*)  have  been  computed  for  known 
m^  and  m,  and  are  presented  in  Table  1.  An  algorithm  for 
using  these  tabled  values  to  find  the  estimates  m^*  and  m2* 
is  presented  in  section  V.B. 

To  describe  t.he  relationship  between  t*,  m^^  and  m.,  ,  a 
three-dimensional  plot  was  made  (Fig.  1)  for  m^^  and  m2  in  the 


range  . 

OK. 01). 09,  . 

l(.l) .9, 

1(1)9.  Fig. 

2  IS  a  corresponding 

plot  of 

EFF(t*)  for 

the  same 

set  of  values 

for  m^  and  m2. 

Figures 

1  and  2  are 

represen 

rations  of  the 

values  found  in 

Table  1. 

For  given  values  of  m^^  and  m., ,  t*  is  the  value  of  t  which 

m.aximizes  the  efficiency  fu.ncticn  EFF(t).  The  values  of  t* 

are  not  necessarily  continuous  as  m^  and  m^  are  varied.  The 

efficiency  function  determined  by  equation  (4.13)  is  bi.mcdai 

for  certai.n  values  of  m,  and  m.,  and  these  modes  occur  near 

1  ^ 

t.he  endpoints  of  the  interval  [0,1).  The  global  .maximum  of 
EFF(t)  on  this  interval  may  occur  at  either  mode  for  different 
values  of  m^^  and  m2,  which  is  the  cause  for  the  discontinui¬ 
ties  of  t*  as  seen  in  Fig.  1.  The  jump  in  the  surface  occurs 
when  m^^  is  about  3.0  for  all  values  of  m2  >  2.0.  Further,  it 
IS  i.nteresting  to  note  that  t*  is  approximately  zero  for 

<  1.0  and  all  values  of  .m2,  except  for  a  singular  point  at 


13 


=  .09  and  =  .01.  EFF{t*)  appears  to  be  a  continuous 
function  of  and  ^2  as  shown  in  Fig.  2  and  is  very  close 

to  1.0  over  much  of  the  range  for  and  m2  described  above. 

The  method  used  for  determining  t*  is  discussed  in 
section  V.A. 
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V.  CALCULATING  TECHNIQUES 


Estimation  of  the  parameters  m^^  and  m2  of  the  Neyman  Type 
A  distribution  using  the  method  outlined  in  section  IV  in¬ 
volves  solving  a  nonlinear  programming  problem;  specifically, 
finding  a  solution  (m^*,  m2*,  t*)  that  maximizes  equation 
(4,13)  subject  to  the  constraint  equations  (4.1)  and  (4.2). 
Since  this  requires  a  large  a.mount  of  computing  power,  an 
approximating  sequential  method  is  proposed  that  can  be  per¬ 
formed  on  a  calculator.  In  order  to  do  this  the  values  of  t* 
must  be  know'n  for  any  m,  and  m^  pair.  These  t*  values  are 
presented  in  Table  1  along  with  the  associated  maxi.mum  EFF 
for  a  wide  range  of  values  of  m^^  a.nd  m., .  Bivariate  linear 
i.nterpolation  as  discussed  in  Abramowitz  and  Stegun  [8]  can 
be  used  to  approximate  t*  between  the  values  of  m^^  and  m2 
given  in  Table  1. 

A.  METHOD  OF  CO.MPUTING  t*  FOR  TABLE  1 

Table  1  was  prepared  by  fixing  m^^  and  m2  at  constant 
values  in  the  range  .01  to  9.0  and  then  conducting  a  single 
variable  search  on  the  function  EFF  (as  a  function  of  t  only) , 
A  Golden  section  search  procedure  was  employed  that  terminated 
at  the  value  of  t*  equal  to  the  midpoint  of  the  searched  in¬ 
terval  on  which  the  value  of  EFF(t)  was  increasing  but  did  net 
change  by  more  than  10  ^  across  the  interval.  Because  of  the 
bimcdality  of  EFF(t)  on  some  regions,  a  search  was  initiated 
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from  the  right  end  of  the  interval  [0,1)  and  a  separate  search 
started  from  the  left  end.  The  global  maximum  was  then  found 
by  comparing  the  resulting  local  maximums. 

3.  METHOD  OF  SOLVING  FOR  m^^*  AND  m2* 

Given  a  sample  of  n  observations  on  the  random  variable  X 
assumed  to  have  an  underlying  Neyman  Type  A  distribution  of 
which  x^  represents  the  ith  sample  observation,  the  algorithm 
for  solving  equations  (4.1)  and  (4.2)  to  find  the  estimators 
m^*  and  m2*  of  the  parameters  m,  and  m2  consists  of  t.he  fol¬ 
lowing  steps: 

STEP  1  -  Compute  the  Method  of  Moments  estimators  from 

the  following  equations 

_  n 

X  =  I  Z  X , 
n  1  =  1 

->  ^  _  -I 

S‘  =  :  (X  -  x)‘  (5.1) 

n-1  1=1 


X 


STEP  2  -  Find  t*  in  Table  1  corresponding  to  mj^  and  m2. 

Note:  If  t*  =  0.  and  the  sample  contains  no  zero  values,  set 

mj^*  =  m^^  and  m2*  =  ^>2  then  STOP. 

STEP  3  -  Solve  the  system  of  equations  (4.1)  and  (4.2) 

for  m, •  and  m.*  with  t  =  t*  by  using  a  Newton-Raphson  pro- 
cedure.  Solving  equation  (4.1)  for  mj^  and  substituting  this 
into  equation  (4.2)  yields 
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_x(g-l)  -  InCt'’^)  =  0  (5.2) 

^2 

where 

m, (t-1) 
g  =  e  “ 

Applying  the  Newton-Raphson  met.hod  on  equation  (5.2)  produces 
as  the  ()c+l)st  iterate  of  1^2* 


(  1 


1  -  e  .t.  1  n  ( t  ) 

jc 


m.  (t-i; 


1  *  e 


■  •'^k  ■ 


where  for  ease  of  typing  nt  =  n.,  and  n.,  =  m_  (the  Method 

-(k)  ^(1)  ^ 

of  Moments  estimator)  .  Whe.n  termination  of  the  above  iteration 

sc.hem.e  occurs  by  t.he  appropriate  criteria  as  set  by  the  user, 

set  m^*  =  m^  and  compute  m, *  from  equation  (5.1). 

^(n)  ^ 

STEP  4  -  If  m-,*  -  m2  <  c  STOP,  where  £  is  set  by  the 

user.  Otherwise  return  to  STEP  2  after  replacing  ^2  with  m2*. 

The  convergence  of  t.his  algorithm  depends  on  t.he  sample 
sice  n  and  how  well  t.he  Method  of  Moments  estimators  approxi¬ 
mate  t.he  true  values  of  m,  a.nd  m_ .  Experience  has  shown  t.hat 
cycling  between  estimates  can  occur  when  the  .11^^*  and  m2* 
values  are  close  to  the  points  where  t*  is  discontinuous  as 
determined  from  Fig.  1.  This  behavior  was  observed  when  the 
above  algorithm  was  applied  to  data  of  sample  sices  as  large 
as  50.  For  most  cases  t.he  estimates  obtained  by  applying  the 
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r 


algorithm  converge  after  two  or  three  iterations.  It  should 


1 


be  noted  from  STEP  2  that  if  t*  =  0  and  the  sample  contains 
no  zero  values,  then  the  .Method  of  .Moments  esti.mators  are  the 
default  values  obtained  when  using  the  algorit.hm. 

C.  .METHOD  OF  DETERMINING  0 

The  value  of  Q  as  defined  in  equation  (3.4)  was  computed 
according  to  equation  (3.5)  for  the  purpose  of  determining  the 
asymptotic  efficiency  of  the  EPGF  esti.mators  given  *n  Table  1. 
The  -xth  power  moment  of  a  Poisson  distributed  random  variable 
(denoted  by  in  equation  ^3.5))  was  computed  by  finding  t.he 


.maximum 

term  of  the 

series  and  using  this 

as  a  scaling  factor 

" 

to  keep 

the  partial 

sums 

wit.hi.n  t.he  li.mit 

of  an  :3.M  360/67 

computer 

,  which  was 

used 

for  all  ccmputat 

ions.  Actual  com.pu- 

of  EFFtt*)  was  computed  as  a  value  greater  than  1.0,  which  is 
i.Tipossible.  It  is  felt  by  this  author  that  in  these  occur¬ 
rences  the  value  of  Q  was  not  computed  accurately  enough. 

The  reason  for  this  is  that  when  m,  or  m-,  was  large  over  300 

terms  were  required  in  t.he  partial  sum  of  Q  in  -der  to  reac.h 

the  stopping  rule  criterion  given  above.  However,  since  the 
size  and  frequency  of  t.hese  inconsistencies  were  small,  it  is 
felt  that  the  other  results  presented  in  Table  1  are  not  com¬ 
promised.  This  IS  justified  by  noting  the  continuous  behavior 
of  EFF(t)  as  given  by  equation  (4.13)  and  observi.ng  the  few 
values  of  EFF(t*)  >  1.0  i.n  Fig.  2.  It  is  net  known  why  these 

occurrences  did  not  happen  at  the  e.xtreme  values  of  m^  =  9.0 

and  =  .01.  As  shown  in  Fig.  2,  t.hey  appear  tc  be  isolated 
points  of  irregularity  in  the  computation  of  EFF{t). 


VI .  COMPARISON  WITH  ESTIMATION  BY  THE  METHOD  OF  MOMENTS 

To  determine  the  utility  of  the  EPGF  method  of  estimation 

for  the  Neyman  Type  A  distribution,  a  comparison  with  Method 

of  Moments  estim.aticn  was  made.  This  comparison  was  performed 

by  computin'^  the  ratio  of  the  asymptotic  efficiencies  of  the 

two  schemes.  The  estimators  of  m,  and  m,  bv  the  Met.hod  of 

1  2  - 

.Mom.e.nts  are  fou.nd  by  solving  the  equations  cf  (5.2)  and  is 
actually  done  i.n  STEP  1  of  t.he  seque.ntial  algorithm,  presented 
in  section  V.3.  The  EPGF  method  is  a.n  attem.pt  to  improve  on 
t.he  Method  of  Mom.ents  estim.ation  of  m.,  and  m,,.  The  amount  of 
im.prove.m.ent  is  determined  by  t.he  increase  in  asym.ptotic  ef¬ 
ficiency  gained  by  applying  t.he  EPGF  method  relative  to  that 
of  the  Method  of  Moments.  In  order  to  find  this  increase. 


m.ptotic  effic 

lency 

of  t.he  Met.hod 

or  the  Neym.an 

Type 

A  distribution 

this  was  done 

and 

IS  shown  below. 

determ.inants 

of  A 

and  C  for  the 

determ.ined  in  a  manner  si.milar  to  section  IV.  3y  using  t.he 

the  following  results  can  be  obtained: 

All  »  -m2 

A, ^  =  -m, 


^21  '^2  " 


■m,  -  2m,. m, 
1  1  2 


equations  in  (5.2) 


thus 


;a1  =  . 


(6.1) 


Also 


*  VAR(.X^)  -  .11^^.712(1  .112 ) 

C,,  =  E(X.^)  -  3E(X,^)E(XJ  ♦  2  (  E(X,1 

X  ^  1  1  J  I 


*  3.11.,  *  ^2  ) 


*“21  ~  1 


Fron  formulas  Jeriv'ed  i.".  (5], 


*  E{X.  -  X)  -  I  V..\R{X,)  ] 


=  11,11^  (  211,1',-,  (1  *  11-,)“  -  (1  -  "n.,  -  6ii- 


30  that 


ii,“n.,“  ;  2.ii,.ii-,vl  T.-,)^  -  -T.-,  211^“  1  n-,"*  j.  (6.2) 


Suostituting  equations  (6.1)  and  (6.2)  into  equation  (3.1) 
yields 


(6.3) 


,  \  [  211,11-,  (1  *  n,) 


where  EFF_^  represents  asynptotic  efficiency  of  the  .Method  of 
rr^ir. 

.Mcnents  esti.T.ators  n,  and  ii-, .  The  ratio  LFF/EFF  fomed  from 

i  2  inn 

equations  (4.13)  and  (6.3)  represents  the  relative  merit  of 
the  EPGF  to  .Method  of  .Moments  es 1 1. ma 1 1 cn  . 
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Values  of  t*  and  EFF/EFF  are  tabulated  in  Table  2  for 

nun 

che  same  range  as  covered  by  Table  1.  This  comparison  of 
the  two  methods  is  most  readily  observed  in  Fig.  3,  which  is 
a  three-di.mensional  representation  of  Table  2.  The  EPGF 
met.hod  is  better  than  the  Method  of  Moments  in  producing 
estimators  of  high  asymptotic  efficiency  for  the  range  of 
parameter  values  covered  in  Table  2.  This  is  clearly  demon¬ 
strated  by  Fig.  3.  The  EPGF  met.hod  reduces  to  the  method  of 
estimation  by  zero  frequency  when  t*  =  0.  In  t.his  case  the 
contours  of  Fig.  3  can  be  fou.nd  in  [4]  and  [9]  . 


VII 


NUMERICAL  EXAMPLES 


EXAMPLE  1 

To  illustrate  the  EPGF  method  and  in  particular  the  al¬ 
gorithm  of  section  V.B,  the  data  of  insurance  claims  caused 
by  accident  or  sickness  cited  by  Shenton  [3]  are  used  and 
reproduced  below.  The  data  are  assumed  to  be  observations 
from  a  Neym.an  Type  A  distribution  and  three  m.ethods  of  esti¬ 
mating  the  param.eters  m,  and  m.,  are  considered; 

1)  the  Method  of  Moments, 

2)  t.he  EPGF  method,  and 

3)  the  Method  of  .Maxim’um  Likelihood. 

The  data  are: 


Value 

0 

1 

3 

•i 

5 

6 

a 

9 

10 

11 

12 

13 

14 

15 

sample 

sample 

sample 


r  requency 

18~ 

185 
200 
16  4 
107 
68 
49 
39 
21 
12 
11 
2 
5 
2 
3 
1 

Size  *  1056 

mean  =  2.8059 
variance  »  6.4106. 
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When  the  EPGF  algorithm  of  section  V.B  is  applied  to  this 
data,  the  following  results  are  obtained; 

a)  from  STEP  1  the  Method  of  Moments  estimators  are 
computed  to  be 

=  2.184  and  m2  =  1.285 

b)  from  STEP  2  the  t*  value  corresponding  to  m^  a.nd  m2 
is  found  by  bivariate  linear  interpolation  in  Table  1  as 
follows ; 

for  m^^  =  2.0  and  m-,  =  1.  285  linear  interpolation  gives 
t*  =  .0'48, 

for  m,  =  3.0  and  m2  =  1.285  linear  interpolation  gives 
t*  =  . 3774 ,  and 

for  =  2.134  and  m.,  =*1.235  linear  interpolation  gives 
f  =  .  1305 

o)  from  STEP  3  t.he  solution  of  (5.2)  yields 
m,  •  =  2.638  and  m.,*  =  1 .064 

where  termination  of  t.he  .Newton- Raphson  met.hod  was  set  for  an 
absolute  error  of  .0001 

d)  from  STEP  4,  where  .  =•  .01,  a  second  application  cf 
the  algorit.hm  is  recuired  and  the  above  m,  •  and  m-,*  is  used 
as  the  new  starting  point  in  STEP  1.  This  second  iteration 
then  yields 

t*  »  .2736  ,  m,*  »  2.609,  and  m,*  =«  l.O^^e. 

Since  the  termination  orite*ia  of  STEP  4  is  not  met  a  third 
Iteration  is  performed  resulting  in  the  following 
t*  *  .2653,  m^^*  a  2.611,  and  m2*  = 
wit.h  the  termination  criteria  being  satisfied 


1.074 


The  computation  of  Maximum  Likelihood  estimators  for  the 
Neyman  Type  A  distribution  is  very  complicated  and  only  ap¬ 


proximate  Maximum  Likelihood  estimates  for  the  above  data 
are  given  in  13).  To  get  an  idea  of  the  computational  effort 
involved,  these  estimators  are  only  the  first  iteration  re¬ 
sults  and  required  several  hours  on  a  desk  calculator  in 
]  949. 

Maximum  Likelihood  estimators  for  the  above  data  were 
found  using  a  TI-59  programmable  calculator  and  the  method 
presented  in  [6].  The  range  of  observed  values  was  small 
enough  so  that  data  storage  was  not  exceeded  on  the  calculator. 
T.he  following  results  were  obtained  in  approxim.ately  t.hirty 
minutes  of  com.puting  time. 

Iteration  ^1 

1  2.475  1.134 

2  2.515  1.116 

3  2.516  1.115 

4  2.516  1.115. 

In  comparison,  the  computing  ti.me  used  to  perform  the 
EPGF  m.ethod  was  approximately  five  minutes  and  the  data 
storage  requirement  consisted  of  only  keeping  eight  numbers, 
which  were  not  influenced  by  the  data  values. 

EXAMPLE  2 

Computer  simulation  was  used  to  generate  eighty-one  data 
sets  of  random  numbers  from  the  Neyman  Type  A  distribution. 
Computer  program  4  in  Appendix  A  was  written  to  simulate  the 
model  presented  in  section  II  and  the  data  sets  produced  are 
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given  in  Appendix  B.  This  program  uses  a  random  number  gen¬ 
erator  that  provides  random  deviates  from  a  Poisson  distri¬ 
bution  and  is  a  local  subroutine.  For  each  data  set  produced, 
the  sample  size  was  fifty  and  the  range  of  values  used  for  .11^^ 
and  m^  was  1(1)9. 

The  EPGF  esti.mates  for  each  data  set  were  found  by  apply¬ 
ing  the  algorith.m  presented  in  section  V.B.  The  t*  value 
required  i.n  STEP  2  was  computed  by  actually  maxi.mizing  the 
efficiency  equation  (4.13)  with  t.he  estimated  values  of  m^ 
and  m, .  The  results  of  t.hese  computations  are  found  in  Table 
3.  The  termination  criteria  of  STEP  4  was  set  at  c  =  .0001, 
so  that  more  iterations  were  required  than  would  normally  be 
expected . 

In  only  two  of  the  eighty-cne  cases  considered,  no  im¬ 
provement  was  achieved  over  the  .Met.hod  of  .Moments  esti.mators. 
These  two  occurrences  are  noted  in  Table  3.  For  the  case 
m,  =  4.0  and  m^  =  T.O,  the  special  term.inaticn  of  STEP  2  of 
the  algorit.hm  was  invoiced;  i.e.,  t*  =  0  and  t.he  sample  con¬ 
tained  no  zero  values.  For  the  case  m,  =  5.0  and  m-,  =  6.0, 

i.  ^ 

cycling  between  t.ne  following  esti.mate  values  occurred. 


.0000  3.2191 

.9310  2.9465 


The  successive  esti.mate  pairs 
of  t*  which  are  on  either  side  of 
seen  in  Fig.  1  of  section  IV.  As 
efficiency  function  is  bi.modal  in 


m. 


9.6798 

10.5754. 

m^^*  and  m^*  produce  values 
the  discontinuous  ridge 
previously  discussed,  the 
this  region  and  causes  the 
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above  behavior.  Since  the  asymptotic  efficiency  is  high  in 
this  region,  a  larger  sample  size  producing  better  initial 
estimates  should  resolve  the  inconsistency  which  results 
when  cycling  occurs. 
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VIII  . 


TABLES 


Tables  1,  2  and  3  mentioned  in  the  previous  sections  are 
contained  on  the  following  pages.  Table  1  presents  the  op¬ 
timal  t  and  efficiency  values  found  for  the  ra.nge  of  m^  a.nd 
m^  from  .01  -  .09  (increments  of  .01),  .1  -  .9  (increments 
of  .1),  and  1.0  -  9.0  (increments  of  1.0).  Table  2  presents 
t.he  optimal  t  and  efficiency  (relative  to  t.he  efficiency  of 
the  Method  of  Moments)  values  found  for  the  sane  range  as 
used  in  Table  1.  Table  3  prese.nts  the  results  of  applyi.ng 
the  algorithm  of  section  V.3  to  t.he  data  sets  of  Appendix  B. 
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IX.  CONCLUSIONS 


The  EPGF  method  provides  estimators  with  asymptotic  ef¬ 
ficiency  close  to  that  of  the  .Maximum  Likelihood  estimators 
over  most  of  the  para.meter  values  considered.  The  sm.allest 
value  found  was  54^  (at  m^^  =  3.0  and  m2  =  9.0)  and  in  m.ost 
instances  the  asymptotic  efficiency  is  greater  than  97» . 

The  EPGF  method  outperforms  the  Method  of  Moments  in  pro¬ 
ducing  estimators  of  consistently  higher  asymptotic  efficiency 
for  all  param.eter  values  from  .01  to  9.0. 

The  efficiency  of  the  EPGF  method  relative  to  the  Method 
of  .Mom.ents  is  less  than  110*  when  m..,  <  .3  for  most  values  of 
m,  .  As  m.ay  be  seen  from.  Table  3,  the  relative  efficienc'. 
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Method  of  Zero  Frequency.  The  EPGF  method  becomes  an  exten¬ 
sion  of  the  Zero  Frequency  .method  when  the  efficiency  of  the 
estim.ators  of  t.he  latter  decreases.  This  happens  when  t*  > 

0.  which  IS  always  the  case  for  m^  >  3.0. 

Although  the  EPGF  method  is  always  better  than  the  .Method 
of  Zero  Frequency  in  terms  of  efficiency,  it  involves  a 
greater  computational  effort  when  t*  *  0.  As  shown  in  [4], 
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the  Method  of  Moments  estimators  have  higher  asymptotic  ef¬ 
ficiency  than  those  of  the  Method  of  Zero  Frequency  whenever 
m^^  >  3.5  for  all  values  of  m^ .  The  efficiency  of  the  .Method 
of  Moments  estim.ators  is  over  20%  higher  than  the  Zero  Fre¬ 
quency  esti.mators  when  m^  ^  4.5.  Combined  wit.h  the  above 
discussion  of  the  EPGF  .method,  whenever  t*  >  0.  the  EPGF 
method  will  yield  estimators  of  substantially  higher  asymp¬ 
totic  efficiency  than  those  found  oy  the  Zero  Frequency  .met.hod. 
The  choice  of  which  method  to  use  always  rests  with  the  user. 
However,  if  it  is  assumed  that  m2  >  . 3  and  that  both  m,  and 
m2  are  not  large  (greater  than  about  7.0)  then  the  extra  ef¬ 
fort  required  by  the  EPGF  m.et.hod  is  compensated  for  by  a 
better  than  10%  increase  in  asymptotic  efficie.ncy  of  the 
estimators.  When  it  is  supposed  that  m^  and  m...  are  not  to  be 
fou.nd  in  the  above  region  t.hen  esti.mation  by  t.he  Method  of 
Moments  is  recomme.nded ,  as  little  increase  i.n  efficiency  can 
be  expected  by  using  t.he  EPGF  .method;  i.e.,  gci.ng  beyond 
STEP  1  of  the  algorit.hm  presented  in  section  V.B. 

As  noted  in  section  VII,  the  EPGF  method  as  applied  by 
the  algorit.hm  of  section  V.B  can  fail  to  converge.  This  is 
characterized  by  cycli.ng  between  two  different  pairs  of 
esti.mators  when  the  sample  size  is  small.  when  this  occurs 
t.he  remedy  is  to  use  a  larger  sample.  The  variance  of  a 
population  having  the  Neyman  Type  A  as  an  underlying  distri¬ 
bution  grows  large  as  m^^  and  m2  increase,  especially  when 
both  are  greater  than  1.0.  Assuming  t.hat  t.he  sample  size  is 


large  enough  for  the  Central  Limit  theorem  to  provide  a  good 
approximation  by  a  Normal  distribution,  the  sample  size  re¬ 
quired  to  provide  a  reasonable  confidence  interval  containing 
the  true  mean  is  quite  large.  Since  x  =  m^m^  is  one  of  the 
estimating  equations  used  in  the  EPGF  method,  it  may  be  dif¬ 
ficult  to  provide  good  initial  estimates  for  the  iterative 
procedure  when  n  is  small.  It  has  been  found  by  e.xperience, 
as  the  sample  size  is  increased,  the  EPGF  esti.mates  tend  to 
cluster  m.ore  tightly  than  the  .Method  of  Mom.ents  estimators 
about  the  true  parameter  values. 

The  EPGF  method  is  not  as  easy  to  use  as  the  .Method  of 
.Moments,  but  wit.h  the  aid  of  Table  2  and  t.he  use  of  the 
algcrit.hm  presented  it  can  be  applied  readily  with  current 
programmable  calculators.  This  is  not  always  the  case  for 
estim.ation  by  .Maximu.m.  Lixeiihccd. 

.Alt.hough  the  EPGF  .method  appears  tc  be  ti.me  consuming 
and  tedious,  it  is  far  easier  to  implement  than  the  Maxi.mum 
Lixeli.hocd  m.ethod.  Dy  using  Table  1,  reasonably  good  esti¬ 
mators  can  be  obtained  with  less  effort  than  required  by  t.he 
latter  method. 
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APPENDIX  A 


COMPUTER  PROGRAMS 

Four  computer  programs  were  used  ir.  the  compilation  of 
this  paper  and  are  given  on  the  following  pages.  All  were 
programmed  in  the  standard  Fortran  IV  language  and  in  double 
precision.  Program  1  was  used  to  compute  the  values  of  the 
i.nformation  determinant  A '  of  section  III  for  m,  and  m_,  in 
the  range  of  Table  1.  Program  2  was  run  to  compute  the 
values  of  t*  and  EFF{t*)  and  produce  Table  1.  Table  2  was 
produced  by  slightly  modifying  program  2  to  compute  and  table 
the  efficiencv  ratio  EFF/EFF  cf  section  VI.  Program  3  was 
used  to  perform  the  algorithm  of  section  V.3  and  produce  the 
esti.T.ates  for  the  numerical  examples  cf  section  VII,  which 
are  found  in  Table  3.  Program.  4  was  used  i.n  si.mulating  ran¬ 
dom  numJsers  from  t.he  Neyman  Type  A  distribution.  The  sub¬ 
routine  LPOISl  found  in  program.  4  is  a  local  random,  number 
generator  for  Poisson  distributed  random  deviates. 

Programs  1,  2,  and  3  are  complete  and  may  be  used  di¬ 
rectly.  Program  4  must  be  modified  to  be  com.patible  with 
the  user’s  local  random  number  generation  subroutines. 
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( i  n  o  o  o  o 


COMPUTER  PROGRAM  1 


C  COMPUTATION  OF  INFORMATION  DETERMINANT  GIVEN  Ml  AND  M2 

IMPLICIT  REAL*8 (A-H,0-Z) 

REAL *8  Ml, M2 
COMMON  Ml , M2 , BETA , DET 
DIMENSION  PARM(27) ,DETQ(27, 27) 

DATA  PARM/l.D-2 , 2 . D-2 , 3.D-2 , 4 .D-2 , 5. D-2 , 6 .D-2 , 7 .D-2 , 

*  8  .  D-  2 , 9  .  D-  2  ,  1 .  D- 1 , 2  .  D- 1 , 3  .  D-  1 , 4  .  D-  1 , 5  .  D- 1 , 6  .  D-  1 , 7  .  D- 1 , 

*  8  .  D-  1 , 9  .  D-  1 ,  1 .  DO  ,  2  .  DO  ,  3  .  DO  ,  4  .  DO  ,  5  .  DO  ,  6  .  DO  ,  7  .  DO  ,  8  .  DO  , 

*  9. DO/ 

CREATE  A  TABLE  OF  DETERMINANT  VALUES  FOR  VARYING 
PARAMETERS  CF 
Ml  AND  M2 
J1  =  27 
J2  =  2  7 

DO  999  I  =  27,27 
Ml  =  PARM(I) 

DO  99  J  =  J1,J2 
M2  =  PARM(J) 

BETA  =  M1*DEXP(-M2) 

VALUES  HAVE  BEEN  SET  FOR  Ml  AND  M2 

COMPUTE  INFORMATION  DETER.MINANT 
CALL  DETCS 
DETQ{I,J)  =  DET 
99  CONTINUE 

IF  ( I .NE. 1)  GO  TO  150 
WRITE  (6,101)  (PARM(L) ,L=J1, J2) 

101  FORMAT(5X,3(3;<,F12.10)  ) 

150  WRITE  (6,151)  Ml ,  ( DETO ( I , D  , L= J 1 , J2 ) 

151  FORMAT(F5.2, 3(3X,F12.10) ) 

999  CONTINUE 

STOP 

END 

C 

SUBROUTINE  DETCS 
IMPLICIT  REAL*8 (A-H, D-Z) 

REAL *8  Ml, M2 
COMMON  Ml, M2 , BETA, DET 
C  COMPUTE  Q  GIVEN  Ml  AND  M2 

H  =  l,D-7 

PO  =  DEXP(3ETA  -  Ml) 

RX  =  BETA 
CX  =•  BETA*RX 
Q  =  QX 

ELN  =  DLOG(BETA)  -t-  BETA 
DO  10  IX  =  1,500 
X  *  DFLOAT(IX) 

ELNEXT  »  SERIES (X+1. DO) 
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RNEXT  =  ESXP(ELNEXT  -  ELN) 

ELN  =  ELNEXT 

QX  =  QX*M2* (RNEXT/X) * (RNEXT/RX) 

RX  =  RNEXT 

IF(QX/Q  .LT.  H)  GO  TO  20 
K  =  IX 

10  Q  =  Q  +  QX 

20  Q  =  M2*M2*PO*Q 

A  =  M1*M2*M2*M2 
B  =  Ml*M2*M2*  (M1+M2-^M1*.M2) 

DET  =  (l.DO  +  .M2)*Q/A  -  B/A 
RETURN 
END 
C 

FUNCTION  SERIES (X) 

IMPLICIT  REAL'S (A-H,0-Z) 

REAL'S  Ml, M2 

COMMON  M 1 , M2 , 3  ET A , D ET 

SERIES  =  O.DO 

C  FIND  INDEX  OF  THE  yjiX  TERM:  AIM 

AIM  =  l.DO 
Cl  =  l.DO 
C2  =  2. DO 
DO  10  K  =  1,100 

FNUM  =  C2'AIM'DLCG  (BETA,  AIM)  C2'X  -  Cl 
FDEN  =  C2'(X  ♦  AIM)  -  Cl 
AINEXT  =  AIM* (Cl  *  FNUM/ FDEN) 

IF(AINEXT  .LE.  1 . D- 5 )  GO  TO  11 
AIM  =  AINEXT 

10  CONTINUE 

11  IF (AIM  .LE.  l.DO)  AIM  =  l.DO 

IF  (AIM  .GT.  l.DO)  AIM  =  DFLOAT  ( I D INT  ( AIM  .  5D0  )  ) 

IM  =  IDINT(AIM) 

C  COMPUTE  .MAX  TERM  OF  THE  SERIES 

ALNIM  =  O.DO 
DO  20  I  =  1,IM 
AI  =  DFLOAT (I) 

20  ALNIM  =  ALNIM  DLCG(AI) 

ALNIM  =  X*DLCG(AIM)  -  AIM'DLOG (BETA)  -  ALNIM 
C  COMPUTE  EACH  TERM  (K)  OF  SERIES  AND  SUM  TO  THE  SERIES 

DO  40  K  =  1,300 
AK  =  DFLOAT (K) 

ALNK  =  O.DO 
DO  30  J  =  1,K 
AJ  =  DFLOAT (J) 

30  ALNK  =  ALNK  *  DLOG(AJ) 

ALNK  =  X*DLOG(AK)  *  AK'DLOG ( BETA)  -  ALNK 
BK  =  ALNK  -  ALNIM 

IF(3K  .GE.  -161. IDO)  SERIES  =  SERIES  DEXP(BK) 

IF(BK  .LT.  -161. IDO  .AND.  K  .GT.  IM)  GO  TO  50 
40  CONTINUE 

C  COMPUTE  THE  NATURAL  LOG  OF  THE  XTH  MO.MENT. 

50  SERIES  =  ALNIM  DLCG  (SERIES) 

RETURN 

END 
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COMPUTER  PROGRAM  2 


C  COMPUTATION  OF  OPTIMAL  T  AND  EFFICIENCY  GIVEN  Ml  AND  M2 

IMPLICIT  REAL*8 (A-H,0-2) 

REAL *8  Ml, M2 

COMMON  Ml,M2,DET,TSTAR,ESTAR 

DIMENSION  PARM(27) ,TOPT(27,27) ,EOPT(27,27) ,DETQ{27,27) 

C  CREATE  A  TABLE  OF  T-OPTIMAL  VALUES  FOR  VARYING 

C  PARAMETERS  OF 

C  Ml  AND  M2 

C 

C  READ  IN  Ml  AND  M2  AND  CORRESPONDING  INFORMATION 

C  DETERMINANT  VALUES  FROM  A  DISK 

C 

J2  =  0 

DO  99  K  =  1,9 

J1  =  J2  1 

J2  =  J1  2 

DO  99  I  =  1,27 

IF (I  .NE.  1)  GO  TO  50 

READ  (5,101) (PARM(J) ,J=J1,J2) 

50  READ  (5,101)  (DETO ( I , J) , J=J1 , J2 ) 

101  FORMAT(5X, 3(3X,F12.10) ) 

99  CONTINUE 

SET  Ml  .AND  M2  VALUES  FROM  PARM 

COMPUTE  OPTIMAiL  T-STAR  =  TOPT  AND  CPTI.MAL  EFFICIENCY 
E-STAR  =  EOPT 
DO  199  I  =  1,27 
Ml  =  PARM (I) 

DO  199  J  =  1,27 
M2  =  PARM(J) 

DET  =  DETQ(I,J) 

CALL  OPT I ML 
TOPT(I,J)  =  TSTAR 
EOPT(I,J)  =  ESTAR 
199  CONTINUE 

C  PRINT  OUT  TABLE  OF  OPTIMAL  T  AND  EFFICIENCY 

lU  =  0 

DO  299  KI  =  1,3 
IL  =  lU  *  1 

lU  =  IL  8 

JU  =  0 

DO  298  KJ  =  1,3 
JL  =  JU  *  1 

JU  =  JL  8 

WRITE  (6,201) 

WRITE  (6,202) 

WRITE  (6,203) 

WRITE  (6,204)  (FARM ( J) , J= JL , JU) 

DO  297  I  =  IL,IU 

WRITE  (6,205)  PARM ( I ) , (TOPT ( I , J ) , J=JL , JU ) 

WRITE  (6,206)  (EOPT  ( I ,  J)  ,  J=rJL,  JU) 

297  CONTINUE 
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WRITE  (6,207) 

WRITE  (6,208) 

WRITE  (6,209) 

298  CONTINUE 

299  CONTINUE 

201  FORMAT( ' 1 ' ,//////) 

202  FORMAT( ’ 0 ', 33X, 'TABLE  1:  OPTIMAL  T  AND  EFFICIENCY 
*  VALUES',//) 

20  3  FORMATCO'  ,9X, 'Ml  VALUES  '  ,  3  5X  ,  '  M2  VALUES') 

204  FORMAT (' O'  , 19X, 9 (2X,F6 . 4) ) 

205  FORMAT ( ' 0 ' , 1 IX , F4 . 2 , 4X , 9 ( 2X , F6 . 4 ) ) 

206  FORMATC  '  ,  19X  ,  9  (2X,  F6 . 4  )  ) 

207  FORMAT (' O' ,/, 22X, ' FOR  EACH  ENTRY  IN  THE  TABLE:') 

208  FORMATCO'  ,21X, 'T*  =  TOP  VALUE  OF  EACH  PAIR') 

209  FORMAT (' 0 ', 21X, ' EFF  =  BOTTOM  VALUE  OF  EACH  PAIR') 
STOP 

END 

C 

SUBROUTINE  OPTIML 
IMPLICIT  REAL*9 {A-H,0-2) 

REAL'S  Ml, M2 

COMMON  Ml , M2 , DET , TSTAR , ESTAR 
FLAG  =  0 

THETA  =  9999. D-4 
E2  =  O.DO 

C  FIND  A  CCARSE  BRACKET  AROUND  THETA:  (A,B). 

5  THETA  =»  THETA  -  l.D-2 

IF (THETA  .LT.  O.DO)  GO  TO  15 
ET  =  EFF (THETA) 

IF  (ET  .LT.  E2)  GO  TO  15 
E2  =  ET 
GO  TO  5 

15  A  =  T!?ETA 

IF  (A  .LT.  O.DO)  A  =  O.DO 
3  =  THETA  -  2.D-2 
:F(3  .GT.  9999. D-4)  3  =  9999. D-4 
C 

PERFCR-M  A  GOLDEN  SECTION  SEARCH  FCR  TSTAR  TO  MJ^X 
C  EFF (TSTAR) 

C  WITH  STARTING  INTERVAL  (A, 3) . 

R  =  (-1.D0  +  DSGRT(5.D0) )/2.D0 
20  A1  =  (l.DO  -  R) • (B  -  A)  ♦  A 
31  =  R*  (B  -  A)  A 
EAl  =  EFF(Al) 

EBl  =  EFF(Bl) 

IF(DABS(EA1  -  EBl)  . LT .  l.D-7)  GO  TO  30 
IF(EA1  .GT.  EBl)  3  =  31 
IF(EA1  .LT.  EBl)  A  =  Al 
GO  TO  20 

30  TSTAR  =«  (Al  31)/'2.D0 
ESTAR  =  EBl 

CHECK  FOR  A  MAX  FROM  THE  LEFT  END 
IF (FLAG  .EQ.  1)  GO  TO  40 
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C 


35 


36 


40 


IF(TSTAR  .LT.  l.D-1)  RETURN 
T  =  TSTAR 
E  =  ESTAR 
THETA  =  O.DO 
E2  =  O.DO 

THETA  =  THETA  +  l.D-2 
IF {THETA  .GT.  99. D- 2)  RETURN 
ET  =  EFF (THETA) 

IF(ET  .LT.  E2)  GO  TO  36 
E2  =  ET 
GO  TO  35 

A  =  THETA  -  2.D-2 

IF (A  .LT.  O.DO)  A  =  O.DO 

3  =  THETA 

FLAG  =  1 

GO  TO  20 

IF(E  .LT.  ESTAR)  RETURN 

TST.^iR  =  T 

ESTAR  =  E 

RETURN 

END 


FUNCTION  EFF(T) 

I.MPLICIT  REAL*9  (A-H,0-Z) 

REAL*  8  .Ml,  M2 

CCM.MON  Ml ,  M2  ,  DET  ,  TSTAR,  ESTAR 

CC.MPUTE  PROBABILITY  GENERATING  FUNCTICN  G(T) 

F  =  DEX?(M2*T  -  .M2) 

G  =  DEX?(-M1* (l.DO  -  F) ) 

COMPUTE  G(T**2) 

F2  =  DEX?(M2*T*T  -  M2) 

G2  =  DEX?(-M1* (l.DO  -  F2)) 

COMPUTE  NUMERAiTOR  OF  EFF(T) 

."N  =  l.DO  *  (M2*T  -  M2  -  l.D0}*F 

FN  =  .Ml*G*G*FN*FN 

COMPUTE  DENOMINATOR  OF  EFF (T) 

FD  =  l.DO  -  M2  -  M1*M2*{T*F  -  l.D0)*(T*F  -  l.DO) 
FD  =  M2* ((l.DO  ^  M2)*G2  -  G*G*FD) 

FD  =  DET*FD 

COMPUTE  EFFICIENCY  EFF(T) 

EFF  =  FN/FD 

RETURN 

END 


68 


o  o 


COMPUTER  PROGRAM  3 


ESTIMATION  OF  Ml  AND  M2 

IMPLICIT  REAL*8 {A-H,0-Z) 

DIMENSION  IX , 10) ,M2 (9) 

COMMON  HATM 1 , HATM2 
READ  (5,9)  (M2 (J) , J=l, 9) 

FORMAT (3X, 915) 

READ  IN  MATRIX  OF  RANDOM  NOS.:  IX 
.READ  (5,10)  (  (IX  { I ,  J)  ,  J=1 ,  10)  ,  1  =  1 , 450) 

FOR.MAT  (13,915) 

NNN  =  50 

.\NN  =  D  FLOAT  I  NNN) 


DO  90  I  =  1,9 
IFLAG  =  IFLAG*(-1) 
IK  =  NNN*  (  :-  1)  -  1 
IKE  =  IK  -  NNN  -  1 


?M1  =  DFLCAT ( IX ( 

IK,  1)  ) 

IF (IFLAG  .LT.  0) 
WRITE  (6,100) 

GO  TC  15 

100 

FOP-MAT  (  '  1'  ,////, 
WRITE  (6,101) 

T  3  7  ,  '  T.A3  LE  3  : 

NUMERICAL  RESULTS') 

101 

FORMATS  O',  ,T19 

■  M  '  '  -'IT  '  V  ' 

,  .  i  A.  0  i  t  ..4. 

,T35,  'Ml'  ,T4  3,  'M.2'  ,T51 

*  ’  T*  •  ,  T5(3,  '  T*  '  ,  T65  ,  'Ml  »  ■  ,  T'i  ,  'M2  *  ’  ,  T90  ,  'NO.  ’  ) 

WRITE  (6,102) 

FOR.MAT  (  '  '  ,T34  ,  '  (M-M)  '  ,  T4  2  ,  '  (.MM)  '  ,T50  ,  '  INIT'  ,T56  ,  'FINAL'  , 

*  T30 , ' ITR. ' ) 

CONTINUE 

DC  90  J  =  2 ,  10 

O  is  —  U  “  ^ 

FM.2  =  DFLOAT  >,M.2  (JK)  ) 

XBAR  =  O.DO 

XSC  =  O.DO 

DC  2  0  K  =  I K  ,  I KE 

XBAR  =  XBAR  -  DFLCAT ( IX ( K , J) ) 

XSC  --  XSC  ♦  DFLCAT  ( IX  (K,  J)  *  IX  (K,J)  ) 

XBAR  =  X3AR/A.NN 

X'.'AR  *  (XSC  -  A.NN*X3AR*XBAR)  /  (ANN  -  l.DO) 

.4ATM2  =  (.XVAR  -  XBAR) /XBAR 

MATMl  =  XBAR/HATM2 

.AMI  =  ?M1 

AM2  =  PM2 

3M1  =  HATMl 

3M2  =  HATM 2 

FIND  ESTIMATORS  USING  EMPIRICAL  FGF 
PAST  =  HATM2 
J  J  J  =  0 
CALL  TOPT  (T) 

IF(JJJ  .EQ.  0)  TI  =  T 

U  vJ  =  <J  J A 
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IF(JJJ  .GT.  20)  GO  TO  65 
TBAR  ^  O.DO 
DO  60  KK  =  IK,  IKE 

IF(IX(KK,J)  .EQ.  0)  TBAR  =  TBAR  1 .  DO 

IF(T  .GT.  l.D-4  .A.ND.  I.X(KK,J)  .NE.  0)  TBAR  =  TBAR  + 

*  T**IX(KK,J) 

CONTINUE 

IF (TBAR  .UE.  l.D-4)  GO  TO  70 
TBAR  =  TB.AR,  A.NN 

CALL  ANEWT  ( TBAR ,  PAST  ,  .XBAR ,  T  ,  PMC  ) 

HATMl  =  XBAR,  PM2 

IF  (DABS  (HATMC  -  P>\2)  .LE.  l.D-4)  GO  TO  80 
HATM2  =  PMC 
GO  TO  30 

WRITE  (6,103)  AMI  ,.A.M2  ,  BMl  ,  3.M2  ,TI 

FCR.MAT  (  ■  0  '  ,  T14,  4F3  .  3,T4  9,F5.  4  ,T56,  ’ - '  ,T64  ,  ' - 

*  T72  ,  ■ - '  ,T9  1 ,  '  -’  ) 

GO  TO  90 

WRITE  (6,104)  AM.1,A.M2,BM1,BM2,TI,T,3.M1,BM.2,JJJ 
FOR.'-IAT  ('0’,T14,4F9.3,T4  9,F5.4,T5  6,F5.4,2F9.3,I5,'*') 
GO  TO  90 

WRITE  i  6  ,  10  5  )  A.M1 ,  .A.M2  ,  BMl ,  3M2  ,  TI  ,  T  ,  HATMl ,  PM2  ,  JJJ 
.=-CR.’‘lAT  l'0',T14,4F9.3,T49,F5.4,T56,F5.4,2F9.3,I5) 
CONTINUE 
WRITE  (6,106) 

FCR.'-IAT  Cl’) 

STOP 

E.VD 


SU3RCUTI);E  ANEWT  ( TBAR  ,  V  IN  ,  XEAR  ,  T  ,  YOUT ) 


IMPLICIT  REAL*9 (A-H,0-2) 

CO.'LMCN  HATM 1 ,  r-ATM2 
TL  »  DLOG(TBAR) 

DO  20  KKK  *  1,20 

GV  =«  DEXP  (YIN*  (T-l.CO)  ) 

FY  =•  iIN*TL  -  X3AR*>1.D0  -  GY) 

FPY  =  TL  -  X3AR*GY*,T  -  l.DO) 

YOUT  »  YIN  -  FY  FPY 

IF (DABS (YOUT  -  YIN)  . LE .  l.D-4)  GO  TO  99 

YIN  =«  YOUT 
20  CONTINUE 

99  RETURN 

END 
C 

SUBROUTINE  TCPT  (T) 

C 

I.MPLICIT  REAL‘9  (A-H,0-Z) 

COMMON  HATM1,HATM2 
FLAG  *  0 

THETA  »  9999. D-4 
E2  =  O.DO 

C  FIND  A  COARSE  BRACKET  AROUND  THETA:  (A, 3) 
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5  THETA  =  THETA  -  l.D-2 

IF (THETA  .LT.  O.DO)  GO  TO  15 
ET  -  EFF (THETA) 

IF{ET  .LT.  E2)  GO  TO  15 
E2  =  ET 
GO  TO  5 

15  A  =  THETA 

IF (A  .LT.  O.DO)  A  =  0 . DO 
3  ^  theta  2.D-2 
IF(3  .GT.  9999.0-4)  3  =  9999. D-4 
R  =  (-1.D0  DSCRT(5.D0)  )/ 2.D0 
20  A1  =  (l.DO  -  R) * (3  -  A)  -  A 
31  =«  R»  (3  -  A)  A 
EAl  =  EFF(Al) 

E31  =  EFF(31) 

IF (DABS (EA1-E31)  . LT .  i.D-15)  GO  TO  30 

:F(EA1  .GT.  E31)  3  =  31 
IF(EA1  .LT.  E31)  A  =  Al 
GO  TO  20 

30  T  *  (Al  -  31)/2.D0 
C  CHECK  FOR  VJ^X  FROM  THE  LEFT  END 

IF (FLAG  . EC •  1)  GO  TO  40 
IF(T  .LT.  l.D-1)  RET’CRi; 

TSTAR  =«  T 
ESTAR  *  E3 1 
THETA  =  O.DO 
E2  *  O.DO 

35  THETA  *  THETA  -  1 . D- 2 

IF (THETA  .GT.  99. D- 2)  RETURN 
ET  *  EFF (THETA) 

IF(ET  .LT.  E2)  GO  TO  36 
E2  =  ET 
GO  TO  35 

36  A  *  THETA  -  2.D-2 

IF (A  .LT.  O.DO)  A  =  O.DO 
B  =  THETA 
F  ..^G  —  I 
GO  TO  20 

40  IF (ESTAR  .LT.EBl)  RETURN 

T  =  TSTAR 
RETURN 
END 
C 

FUNCTION  EFF(C) 

C 

I.MPLICIT  REAL'S  (A-H,0-Z) 

COMMON  HATM1,HATM2 
A.M1  »  HATMl 
A.M2  »  HATM2 
F  =  DEX?(AM2'0  "  A.M2) 

G  =  DEX?(-AM1*(1.D0  -  F) ) 

F2  =  DEXP  (A.M2*'C*Q  -  A.M2 ) 

G2  *  DEXP  (-A.M1' (1  .DO  -  F2)) 


FN  *  L.DO  (AM2*Q  -  AM2  -  l.DO)*F 
FN  =»  AM1*G*G*FN*FN 

FD  =  l.DO  AM2  -*■  AM1*A.M2*  (0*F  -  l.DO)*(Q*F  -  1 .  DO ) 
FD  =»  AM2*((1.D0  -t-  AM2)*G2  -  G*G*FD) 

EFF  =  FN/'FD 

RETURN 

END 


COMPUTER  PROGRAM  4 


C  RANDOM  NUMBER  GENERATION  FROM  NEYMAN ' S  TYPE  A 

C  DISTRIBUTION 

DIMENSION  AD (50) 

DATA  ISEEDl, ISEED2 , I SORT , KS I 2E/4 3 1 2 1 9 , 1 8 76 4 8 , 0 , 50/ 
CALL  OVFLOW 
DO  10  Ml  =-  1,9 
PMl  =  FLOAT (Ml) 

CALL  LPOISl {ISEED1,AD,KSI2E,PM1, I SORT) 

WRITE  (6,1) 

1  FORMAT( ' 1' ,///////) 

WRITE  (6,2) 

2  FORMAT  (■  '  ,  14X,  ■  RA.NDGM  NOS.  FROM  THE  NEY.MAN  TYPE  A 

•  DISTRIBUTION' ) 

WRITE  (6,3) 

3  FORMATCO',  ,',12X'M1  VALUE  '  ,  T  4  0  ,  ’  y.2  VALUES') 

CALL  DGHTER(AD, I SEED2 , Ml , FS I Z E ) 

10  CONTINUE 

STOP 
END 
C 

SUBROUTINE  DGHTER(AD, ISEED2 ,Ml^ 

DIMENSION  AD (50)  ,Y  (100)  , IX (50, 9) 

WRITE  (6,99)  (J,J=1,9) 

99  FOR.'^AvTCO'  ,  19X,9I5,//) 

DC  200  K  1,50 
ISIZE  =•  INT  (AD  (K)  ) 

IFdSIZE  .NE.  0)  GO  TO  53 
DC  50  M2  *  1,9 

50  IX(K,M2)  =  0 

oO  1 0  D 

53  DO  100  M2  =  1,9 

rM2  =  FLOAT (M2) 

CALL  LPOISl ( ISEED2 ,Y, ISIZE, PM2,0) 

IX(K,M2)  =  0 
DO  100  I  =  1, ISIZE 

100  :X(K,M2)  *  IX(K,M2)  INT(Y(:)) 

5  5  WRITE  (6,101)  Ml ,  ( IX  (K , J)  , J= 1 , 9 ) 

101  FOR.MAT(’  '  ,  lOX,  15 , 3X,  915) 

200  CONTINUE 

RETURN 

END 
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